APS /Manuscript 



0(N) Algorithms in TBMD Simulations of the Electronic Structure of Carbon 

Nanotubes 

G. Dereli 

Department of Physics, Middle East Technical University, 06531 Ankara, Turke£\ 

C. Ozdogan 

Department of Computer Engineering, Qankaya University, 06530 Ankara, Turkey^ 

(Dated: February 2, 2008) 

The O(N) and parallelization techniques have been successfully applied in tight-binding MD 
simulations of SWNT's of various chirality. The accuracy of the O(N) description is found to be 
enhanced by the use of basis functions of neighboring atoms (buffer). The importance of buffer size 
in evaluating the simulation time, total energy and force values together with electronic temperature 
has been shown. Finally, through the local density of state results, the metallic and semiconducting 
behavior of (10 x 10) armchair and (17 x 0) zig zag SWNT's, respectively, has been demonstrated. 
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I. INTRODUCTION 

The carbon nanotubes are playing a major role in the 
design of next generation nanoelectronic and nanoelec- 
tromechanical devices due to their novel mechanical and 
electronic properties^. The conductivity behavior of sin- 
gle walled carbon nanotubes (SWNT) is mostly deter- 
mined by the chirality of the tubes. Depending on its 
chirality, a SWNT could be a conductor, semiconductor, 
or insulator. It is now widely known that the conduc- 
tivity of the tubes may also change due to the presence 
of defects as well as radial deformations^. Tight bind- 
ing calculations have shown that deformations such as 
uniaxial compressive/tensile or torsional will also mod- 
ify the band gap of the nanotubes and under such de- 
formations SWNT undergo conducting-semiconducting- 
insulator transitions?' 4 . Real space algorithms has been 
successfully used to perform the ab initio electronic struc- 
ture calculations in the literature^-^. The main objective 
of this paper is the use of O(N) parallel tight binding 
molecular dynamics method in studying the electronic 
structure of SWNT's with diameters going upto 2nm. 
We applied the O(N) and parallelization techniques in 
particular to (10 x 10) and (17 x 0) SWNT's. 

The tight binding theory (TB) has been established 
as a good comprimise between ab initio simulations and 
model-potential ones, bridging the gap between them, 
either as far as the overall numerical efficiency and/or 
as far as the accuracy were concerned. Tight Binding 
Molecular Dynamics (TBMD) is a computational tool de- 
signed to run finite-temperature MD simulations within 
the semi-empirical tight-binding scheme^jiS. This tech- 
nique can be used to simulate material systems at differ- 
ent conditions of temperature, pressure, etc., including 
materials at extreme thermodynamical conditions. The 
electronic structure of the simulated system can be cal- 
culated by a TB Hamiltonian so that the quantum me- 
chanical many-body nature of interatomic forces is nat- 
urally taken into account. The traditional TB solves the 



Schrodinger equation by direct matrix diagonalization, 
which results in cubic scaling with respect to the number 
of atoms ( 0(N 3 )). The 0(N) methods, on the other 
hand, make the approximation that only the local en- 
vironment contributes to the bonding and hence band 
energy of each atom. In this case the run time would 
be in linear scaling with respect to the number of atoms. 
Moreover, O(N) schemes can be efficiently parallelized 
through the use of message passing libraries. The mes- 
sage passing libraries such as PVM and MPI are making 
the simulations possible on clusters of computers. We 
applied these two techniques {O(N) and parallelization) 
succesfully to the SWNT simulations on a cluster of eight 
PC's. Details of the computational study can be found 
in our previous workii^. 



II. THE METHOD 

Traditional TB solves the Schrodinger equation in the 
reciprocal space by direct matrix diagonalization, which 
results in cubic scaling with respect to the number of 
atoms. The O(N) methods on the other hand, solve for 
the band energy in real space and make the approxima- 
tion that only the local environment contributes to the 
bonding, and hence band energy, of each atom. All the 
O(N) methods in which the properties of the whole sys- 
tem are computed (such as the charge distribution, the 
total energy or the forces on all atoms), provide nec- 
essarily approximations to the exact solution of the ef- 
fective one-electron Hamiltonian. These approximations 
are based on physical assumptions, which are generally 
connected to the above mentioned locality or nearsight- 
edness principle in one way or another. Most of the im- 
plementations of the O(N) procedure have been devel- 
oped for the orthogonal tight -binding Hamiltonian. The 
O(N) techniques may be roughly grouped into two cate- 
gories: variational methods and moment-based methods. 
There are two types of variational methods: the density 
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matrix methods and localized orbital methods. There is 
also a variety of moment methods. The O(N) scaling, 
in these approaches, arises from the decay and/or trun- 
cation of these respective quantitie s 13 15 . In our work, 
we adopted a divide and conquer approach (DAC) first 
proposed by Yan g 1 ^ 1 7 as a linear-scaling method used 
to carry out quantum calculations. The basic strategy 
of this method is as follows: divide a large system into 
many subsystems, determine the electron density of each 
subsystem separately, and sum the corresponding con- 
tributions from all subsystems to obtain solely from the 
electron density^. Each subsystem is described by a set 
of local basis functions, instead of the entire set of atomic 
orbitals. The accuracy of the description is enhanced by 
the use of basis functions of neighboring atoms. These 
neighboring atoms are called the buffer. Here the form 
of the Schrodinger's equation of the buffer will be as in 
Refiii. The eigenvalues and vectors are found by di- 
agonalizing the Hamiltonian matrices for each subsys- 
tem. Let M be the number of atoms in the buffer region 
while N the number of atoms in a subsystem and NCell 
the number of subsystems. A subsystem will be labeled 
by a. P % shows the projection of the i'th electron and 
O l = /((£, — /i)/fcsT) the occupation, n will be the 
total number of atoms in the system, MM the number 
of atoms in the buffer region that are in the interaction 
distance (cutoff). Thus we write 
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where H(j, i) is the ji th eigenvector after diagonalization 
scheme and 



O l = 



l + f((£i-v)/k B Ty 



(2) 



Then 



-IN 



l + /((£- M )/fe B T) ^ 
and the subsystem density 

Pa = E^" 



*]T|tf(j,z)] 2 (3) 



The trace 
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trace = p a 



(4) 



(5) 



must be equal to number of electrons in the system so 
that the error 



error = trace — 4 * n. 



(6) 



In the above expressions, f(x) = 1/(1 + exp(xj) is the 
Fermi function, p is the chemical potential for the elec- 
trons, Ub is the Boltzmann constant, and T is the tem- 
perature of the electrons in degrees Kelvin. In case the 



error exceeds the accuracy needed within the desired elec- 
tron temperature, the chemical potential is recalculated 
from 



('n 



-error 



dp 



(7) 



where 

NCell 4JV 

d P = E E * P% ) * C 1 - O l )/k B T] . (8) 
i=i 

This procedure is repeated until the desired level of accu- 
racy is reached. The final value for the chemical potential 
is identified as the Fermi energy level of the system. 
The band structure energy of the system is calculated 

as 

NCell 

Eiy S = ebstot a (9) 

a=l 

where ebstot a shows the contribution of a subsystem to 
the band structure energy of the system: 

4JV / 4AT 

e6stot Q = ^^[| 2* density a {i,j) *TL{i,j) 



+density a (i, i) * 7i(i, i)] 



(10) 



where 



4N I 4N 4W 

density a (k,j) = J} EE ^ * H( - k > * 0% 
j=i yfc=j i=i 

(47V 4Af \ 
J20.5*H(j,i)*H(k,i)*O i \] (11) 
k=4N+l i=l / 

and Ti.(i)j) is the ji th element of the Hamiltonian matrix 
of the subsystem. 

The next step is to find the forces that each atom ex- 
periences arising from the electronic structure, i.e. in the 
x-direction: 

NCell AT 

/*;=,..„= E E •/ ("J 

a—l i—1 

where 

J\[N 4 4 

/-,= E E E density a (4(i -l) 

j — l im—ljm—l 

4(j — 1) + jm) * Force(im, jm) (13) 

and Force(im, jm) has the same form as in Ref^i. Total 
energy of the system has the form, 

Etot — Ebs + E rep (14) 

where E rep has the same form also as in Reft**. The 
energetics and forces are now calculated and then molec- 
ular dynamics scheme is applied and this procedure is 
continued until structural stability is sustained. 
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III. RESULTS AND DISCUSSION 

We applied the O(N) and parallelization techniques 
in particular to (10 x 10) and (17 x 0) SWNT's. Ob- 
viously the results obtained with O(N) algorithm must 
be consistent with 0(N 3 ) results for the same SWNT. 
Buffer size and the electronic temperature are the im- 
portant O(N) parameters that affects the results. We 
produced the Fermi-Dirac distribution, local density of 
states and energetics for the above types of SWNT's. 
Through these, it is possible to distinguish between the 
metallic and semi-conducting behavior of SWNT's de- 
pending on their chiralities. 

The buffer size is an important parameter in the O(N) 
simulations. Each subsystem has its own buffer region so 
that its own small sized Hamiltonian matrix. After di- 
agonalizing this Hamiltonian matrix, the eigenvalues and 
eigenvectors are obtained. The next step is to obtain all 
these informations for all subsystems and then calculate 
the overall system property; chemical potential. This pa- 
rameter gives us the value for Fermi energy level. Then, 
the forces that each atom experiences and the contribu- 
tion of each subsystem to band structure energy of the 
system are calculated. All these procedures are repeated 
through each MD time step. The results obtained with 
O(N) algorithm must be consistent with 0(iV 3 ) results 
for the same system. To ensure this, the value for the 
buffer size parameter must be investigated. 

Another important parameter in the simulation is the 
cuboidal box size. We took the cube size equal to the 
distance between the layers in the tube so that each cube 
has equal amount of atoms. This also provides the same 
number of interacting neighbor atoms (buffer) for each 
subsystem. The PBC is applied in the z-direction only. 
Hence, the system behaves as infinitely long tube. In 
Fig. ^ h is seen that the difference with 0(iV 3 ) total 
energy result for 18 layers and 24 layers are exactly same, 
since PBC works well. We have chosen 20 layers for both 
(10 x 10) and (17 x 0) tube structures for the rest of our 
study. 

Buffer atoms are selected using a distance criterion, Rf,. 
That is, if an atom is within a distance Rb of a subsystem, 
this atom will be included as buffer atom for that subsys- 
tem. The diagonalization for a subsystem is performed 
with atomic basis functions of the subsystem atoms and 
buffer atoms, and the computational effort scales as iV 3 , 
where N a is the number of basis functions in the a sub- 
system and its buffer region. After diagonalization, the 
resulting eigenvalues and eigenvectors give us necessary 
information for local Density of States (1DOS) and for 
force expressions to evaluate the next MD iteration. 

In Fig. ^ the effect of the buffer size on the total 
energy within the given constraints such as boxsize, elec- 
tronic temperature is given. It is seen that the effect of 
the buffer size on the O(N) TBMD is very important. For 
the 10x10 SWNT difference with 0(iV 3 ) TBMD result 
(error) decreases when the buffer size is increased; then 
reaches to desired accuracy and fluctuates around this 



value. Buffer size is important in evaluating the simula- 
tion time, energy and force values. Such as, if the buffer 
size is chosen a higher value than necessary, it will affect 
the simulation time in cubic manner since the Hamilto- 
nian matrix is constructed with respect to the number 
of interacting atoms in the buffer region. On the other 
hand, if this parameter is chosen a low value then it will 
not be able to to produce the correct energy and force 
values. In Figs. [21 and the effect of the buffer size 
on the O(N) TBMD for the (10 x 10) and (17 x 0) tube 
structures together with the effect of the electronic tem- 
perature are also given. From Fig. [3 the buffer size for 
(10 x 10) tube is chosen as 4.9 A and for (17 x 0) tube 
it is 5.7 A ( see Figure [3J. It is important to keep the 
buffer size parameter as small as possible and at the same 
time, it must be able to produce the same values with the 
0(A 3 ) TBMD results. The buffer size for (17 x 0) tube 
converges to desired accuracy much later than (10 x 10) 
tube. This results in much longer simulation time. 

In the Figs. |21 and [31 the effect of the buffer size to- 
gether with the varying electronic temperatures (from 
k B T = 0.000001 eV («0.012 K) to k B T = 0.1 eV 
(«1200 K)) on the O(N) TBMD total energy value for 
the (10 x 10) and (17 x 0) tube structures are also given. 
It is seen that when the buffer size value is small, elec- 
tronic temperature has a slight effect on the energy value. 
These values are the static results without performing 
simulation. The effect of the electronic temperature may 
be impressive during the simulation, when the forces be- 
tween the atoms become dynamic. Therefore during sim- 
ulations, it is safe to choose the electronic temperature 
as room temperature. 

Next the effect of the electronic temperature on the 
total energy is investigated. The energetics for the pris- 
tine (10 x 10) tube with different electronic temperature 
values k B T = 0.025 eV («300 K) and k B T = 0.05 eV 
(«600 K) are studied. Results are given in the Fig. 0] 
The upper graph is for the room temperature and the 
lower is the twice of the room temperature. Having an 
equal average energy value both graphs show similar be- 
havior. They fluctuate around almost the same value. 
This is reasonable because they both simulate the same 
system with different electronic temperatures. The pat- 
tern at the above graph is more dense than the lower one. 
This is because of the hotter electrons in the system. 

In order to observe the effect of electronic temperature 
on the Fermi-Dirac Distribution of the pristine (10 x 10) 
and (17 x 0) tubes with different electronic temperature 
values k B T = 0.01 eV («120 K) and k B T = 0.1 eV 
(wl200 K) are studied. Results are given in the Figs. 
[5] and HJ The upper graphs in the figures are at 120 K 
while the lower ones are at 1200 K. It is observed that 
as the electronic temperature is increased the graphs are 
broadening. Since less electronic state is populated at the 
low electronic temperature condition there is no widening 
for the upper graphs as expected. 

The density of states is obtained from the general for- 
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mula, 

9i 8) = d J§i = ^±Azm (15) 

where N is the number of electrons in the system and 
equal to, 

NCell 4J\f 4N 

(16) 

In the Eq. the statement is that; if there is a change in 
the slope this gives us the information about the existence 
of populated electronic state. The criteria is the change 
in the slope for the Figs. 171 17(71 

In the Figs. [7J-E3 the 1DOS graphs for the pristine 
(lOx 10) and (17x0) tubes for different electronic temper- 
ature values k B T = 0.1 eV («1200 K), k B T = 0.05 eV 
(«600 K) and k B T = 0.05 eV («600 K), k B T = 0.025 eV 
(«300 K); respectively are given. In the Figs. [7Jand[51 
only a selected range for 1DOS are given to better un- 
derstanding of the behavior of electronic states near the 
Fermi-Energy level for the tube structures (10 x 10) and 
(17 x 0), respectively. The other figures, namely |H1 and 
1101 give the same information but in the full range. It 
is seen that when the electronic temperature is increased 
the graphs begin to be smoother since higher amount of 
electronic states are populated. But, the peaks at and 
around the Fermi energy level are at the same positions 
for different electronic temperatures as expected. 

The Fermi energy level values are very similar (around 
3.7 eV) for both tube structures. Although they have dif- 
ferent chirality this is expected because two tubes have 
the same radii. The formula for the DOS gives the elec- 
tronic state population for the different energy values. It 
is found that the (10 x 10) tube has metallic behavior 
since it has states around Fermi energy level and a wide 
band gap but on the other hand the (17 x 0) tube has 
semiconducting behavior since it has no states around 
Fermi energy level and small band gap as expected. 

We have also calculated the band gap values for (10 x 
10) and (17 x 0) tubes as 2.01 eV and 0.53 eV; respec- 
tively. In the literatur e 19 24 , the proposed model values 
are calculated by the formulas 2joa c - c /d and 6joa c ^ c /d 
(where 70 = 2.5 — 2.7 eV, a c _ c = 0.142 nm, and d for di- 
ameter in (nm)) for semiconducting and metallic tubes, 
respectively. A theoretical value for 70 of 2.5 eV has 
been estimated by 2 ^ using a first-principles local density 
approximation (LDA) to calculate the band structure of 
armchair carbon nanotube. As it is discussed inS this 



type of calculations give 10-20% smaller value for arm- 
chair nanotubes. The band gap values for the (10 x 10) 
and (17x0) Carbon Nanotubes using this model are 1.62- 
1.75 eV and 0.54-0.58; respectively. Our O(N) TBMD al- 
gorithm gives good energy band gap result for the (17x0) 
tube. For the (10 x 10) armchair tube the model value is 
lower than our calculated value. On the other hand, the 
behaviors of the local density of states graphs are (see 
Figs. IT HTrj|> as expected. For the 10x10 tube (metallic 
behavior), band gap is wide and there are states popu- 
lated around Fermi energy level and for the (17 x 0) tube 
(semiconducting behavior), band gap is narrow and no 
states around Fermi energy level as expected. 



IV. CONCLUSION 

In this study, the details of O(N) TBMD algorithm 
is given. It is described that how a system is divided 
into many subsystems and how their contributions give 
overall system properties (such as charge density, band 
structure energy) by using nearsightness principle. This 
principle uses the approximation that only the local envi- 
ronment contributes to the bonding of each atom. This 
gives us the opportunity for linear scaling. The main 
problem in the traditional TB is the increasing system 
size. When the system size increases (N), the time to di- 
agonalize the constructed Hamiltonian matrix becomes 
in the order of TV 3 . The O(N) algorithms overcome this 
bottleneck and the behavior has a linear scaling. It is 
shown in Ref. 1 ^ that our O(N) algorithm scales linearly 
for increasing system size. 

To conclude, our methodology is able to produce the 
physical properties such as Fermi-Dirac Distribution, lo- 
cal Density of States and energetics for the Carbon Nan- 
otubes. The structural stability under extreme condi- 
tions such as uniaxial strain will be studied separately. 
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FIG. 1: The difference of 0(iV 3 ) total energy result (- 
8.350497775) with O(N) total energy result for the variation 
of O(N) parameter ( Buffer Size) for the 24 Layers and 18 
Layers 10x10 Tube with box size is equal to 1.229 A , the 
electronic temperature is UbT = 0.005 eV . 




FIG. 2: The effect of Buffer Size for the Total Energy value 
for the tube structure 10x10, with different electronic temper- 
ature values, respectively. 
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FIG. 3: The effect of Buffer Size for the Total Energy value for 
the tube structure 17x0, with different electronic temperature 
values, respectively. 
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FIG. 4: The effect of electronic temperature on the total 
energy (Tube 10x10 no strain T= 300 K) for the values 
k B T = 0.025 eV and k B T = 0.05 eV, respectively. 
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FIG. 5: Fermi-Dirac Distribution Function vs Energy for var- 
ious electronic temperatures for the Tube Structure 10x10. 
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FIG. 6: Fermi-Dirac Distribution Function vs Energy for var- 
ious electronic temperatures for the Tube Structure 17x0. 
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FIG. 7: Local DOS vs Energy for different electronic temper- 
atures for the Tube Structure 10x10.. 




FIG. 8: Local DOS vs Energy for different electronic temper- 
atures for the Tube Structure 10x10 
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FIG. 9: Local DOS vs Energy for different electronic tem- 
peratures for the Tube Structure 17x0. 




FIG. 10: Local DOS vs Energy for different electronic tem- 
peratures for the Tube Structure 17x0. 



